The structure of binary Lennard-Jones clusters: The effects of atomic size ratio 
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We introduce a global optimization approach for binary clusters that for a given cluster size is able 
to directly search for the structure and composition that has the greatest stability. We apply this 
approach to binary Lennard-Jones clusters, where the strength of the interactions between the two 
atom types is the same, but where the atoms have different sizes. We map out how the most stable 
structure depends on the cluster size and the atomic size ratio for clusters with up to 100 atoms 
and up to 30% difference in atom size. A substantial portion of this parameter space is occupied 
by structures that are polytetrahedral, both those that are polyicosahedral and those that involve 
disclination lines. Such structures involve substantial strains for one-component Lennard-Jones 
clusters, but can be stabilized by the different-sized atoms in the binary clusters. These structures 
often have a 'core-shell' geometry, where the larger atoms are on the surface, and the smaller atoms 
are in the core. 

PACS numbers: 61.46.-|-w,36.40.Mr 



I. INTRODUCTION 

There has been much recent interest in binary clusters, 
both from a fundamental and a technological perspective. 
For example, alloy clusters are of particular importance, 
because of their potential catalytic properties ii*S Further- 
more, binary clusters offer the opportunity to tailor their 
properties through the choice of atom types and compo- 
sition, potentially leading to new behaviour that is not 
possible for single component clusters. Here our focus is 
on the structure of binary clusters, which is one of their 
most important properties, and is a prerequisite for un- 
derstanding many of the other properties of such clusters. 

For monoatomic clusters, Lennard-Jones (LJ) clusters 
provide a well-characterized model system, for which the 
effects of a cluster's finite size on the structure, ther- 
modynamic and dynamic properties are well understood. 
For example, putative global minima are now available 
for all LJ clusters with up to 1600 atoms t^i^i^iSiL^ and 
the size evolution of the structure of larger clusters is 
well understoodi^iiSiii 

For binary clusters, there is a similar need for an 
archetypal system to understand the effects that can con- 
trol the structure of a binary cluster. We propose that 
binary Lennard-Jones (BLJ) clusters could provide just 
such a model system, and this paper aims to start the 
systematic exploration of this model. One advantage of 
using BLJ clusters as a model system is the relative sim- 
plicity of the potential. There are only four effective pa- 
rameters that characterize the interactions between the 
different combinations of atoms, making it possible to 
systematically study the structure of BLJ clusters as a 
function of these parameters. By contrast, even the sim- 
plest many-body metal potential would have a consider- 
ably greater number of parameters, and so it only be- 
comes feasible to study a series of example binary metal 
systems, which just represent a set of points in this state 
space of potentials. 



Binary clusters also offer considerable additional 
challenges to the theoretician, compared to the one- 
component case. Firstly, for a given cluster, there are 
many more minima on the potential energy surface, be- 
cause of the presence of "homotops" ^ isomers with the 
same geometric structure, but which differ in the la- 
belling of the atoms. Secondly, the composition provides 
an additional variable to consider. For example, the task 
of obtaining the lowest-energy structures for all compo- 
sitions and all sizes up to 100 atoms would require 5050 
different global minima to be found. Instead, most stud- 
ies have either just considered a few sizes and explored 
how the structure depends on compositionjiii^iiSiiSiii 
or kept the composition fixed and studied the size 
dependence. ^^'^^ By doing such a selective survey, there is 
the possibility that the most interesting and stable struc- 
tures for the systems are missed. 

Here, we take a different approach, neither trying to 
find the global minima for every size and composition, 
nor taking an arbitrary cut through this space. Since 
one is usually just interested in particularly stable struc- 
tures, finding all the global minima is unnecessary, and 
so we instead directly search for these particularly stable 
structures. In particular, we use the composition as a 
variable in our global optimization and for a given size 
we attempt to find the most stable composition. 

Our focus in this paper is to use BLJ clusters to un- 
derstand how different types of cluster structure can be 
stabilized purely through the two atom types being of 
different size. In particular, we are interested in polyte- 
trahedral structures^S for which all the occupied space 
can be divided up into tetrahedra with atoms at their 
corners. For a one-component system such packings are 
said to be 'frustrated', because regular tetrahedra do not 
tesselate, as illustrated in Fig. ^ The best local pack- 
ing of tetrahedra involves packing five tetrahedra around 
a common edge, but with regular tetrahedra there is a 
small angular deficit. As larger structures are made up 




FIG. 1: (Colour online) Examples of the strain involved in 
packing tetrahedra. (a) Five regular tetrahedra around a com- 
mon edge produce a gap of 7.36°. (b) Twenty regular tetra- 
hedra about a common vertex produce gaps equivalent to a 
solid angle of 1.54 steradians. 



of regular tetrahedra, these gaps grow rapidly in size. 
So for a 13-atom icosahedron, which can be thought of 
as a packing of 20 slightly irregular tetrahedra around a 
common vertex, the distance between adjacent atoms on 
the surface is 5.15% longer than that between the cen- 
tral atom and a surface atom. For monoatomic systems, 
there will be an energetic penalty associated with this 
strain. However, the associated strain can be completely 
removed in a binary system simply by choosing the cen- 
tral atom to be 9.79% smaller. Similarly, Frank-Kasper 
phases, bulk polytetrahedral crystals, are only found for 
alloys i^^i^^i^ Furthermore, stabilization of polytetrahe- 
dral structures in binary metal clusters in which there is 
a size mismatch between the two atom types has previ- 
ously been seeniiii^ii^ 

In sectionn]we describe the choice of potential param- 
eters, how we analyse the energetics, and our global opti- 
mization approach. In section llll Al for three cluster sizes 
we present case studies of how the lowest-energy struc- 
tures depend on both the composition and the atomic 
size ratio, and in section llll Bl we present our systematic 
survey of the optimal structures and compositions for all 
clusters with up to 100 atoms and up to 30% difference in 
atom size. A brief report of some of the work presented 
here has appeared previously^ 



II. METHODS 



Potential 



Here, we use a binary Lennard- Jones (BLJ) potential: 
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where a and (3 are the atom types of atoms i and j, and 
tap and i^/^aap are the pair well depth and equilibrium 
pair separation, respectively, for the interaction between 
atoms i and j. In its most general form the BLJ po- 
tential has four effective parameters, namely tAB, ^bb, 
Gab and gbb, if ^aa and gaa are used as the units 
of energy and length, respectively. Here, as we wish to 



consider the effects of purely the size ratio on the most 
stable structures, we choose tAA — ^ab — ^bb — £■ 
Furthermore, as we define gab using the Lorentz rule, 
CAB = {o'aa + crBB)/2, in this case there is effectively 
just one parameter in the potential, namely the size ra- 
tio of the two atoms, ctsb/ctaa. 

Initially, we thought of directly searching for the most 
stable composition and size ratio for a cluster of a partic- 
ularly size in our global optimization runs. However, we 
quickly found that increasing the size disparity between 
the two types of particles leads to a virtual monotonic de- 
crease in the energy. Consequently, our optimization runs 
led to structures with huge differences in the sizes of the 
atoms. These structures consist of a core of tiny atoms 
surrounded by a shell of large atoms, where the large 
atoms are able to interact strongly with all the atoms in 
core. They are clearly unphysical and so this approach 
was abandoned in favour of using only the composition 
as a variable during the global optimization, and consid- 
ering different size ratios independently. 

We chose to look at the structures for size ratios in 
the range 1 < (tbb/<^aa < 1-3, as we wished to see 
if polytetrahedral structures would be stabilized as we 
move away from the one-component Lennard- Jones ref- 
erence system (i.e. ctbb/o'aa = !)■ Due to the sym- 
metry of the energetic interactions, exactly the same 
structural behaviour will be seen in the parameter range 
0.769 < (JBB/crAA < 1, except that the role of the A and 
B atoms will be reversed. Six values were considered; 
namely, cbb/ctaa — 1-05, 1.1, 1.15, 1.2, 1.25 and 1.3. 

There has been a certain amount of previous work 
on binary LJ clusters, as a general model for bi- 



nary systems f 



25.26.27.28.29 



as a model of binary rare-gas 



clustersiSiSSi^iiSSi^ and as a playground for testing how 
the potential parameters could be used to tailor the clus- 
ter's structural properties i22iMi2Si2& Of these studies, a 
few have considered the same set of parameters as we 
have heref2^iSL2S with some tendencies to form core-shell 
clusters notedi2S*2i However, nothing like the systematic 
survey presented here has been tried previously. 

To apply the current approach, one needs to be able 
to compare the relative stabilities of clusters with the 
same number of atoms but different compositions. For 
the current system, as all the eij are the same, the most 
natural way is to compare their absolute energies directly. 

For systems, where the energetic interactions are more 
varied, this issue becomes more subtle. There are a num- 
ber of approaches one might take depending on to what 
the 'stability' is being measured with respect. For ex- 
ample, a simple approach would be to compare the ener- 
gies to a linear interpolation between the energies of the 
pure systems. In bulk, this would amount to compar- 
ing the stability with respect to separation into the pure 
phases. For clusters, it would be sensible to somewhat 
amend this approach, as it would automatically put the 
pure clusters on the same energetic footing, irrespective 
of whether they represented particularly stable sizes (or 
not) for the pure systems. For instance, one could in- 



stead measure the energies with respect to a linear inter- 
polation between smoothly varying functions, call them 
E2vg{N), that captured the general size dependence of 
the energies of the pure clusters. 

Alternatively, one might be interested in maximizing 
the stability of the cluster with respect to the liquid state, 
i.e. finding the composition with the maximum melting 
point. For example, for systems where the cross-terms 
favour mixing, there will be a stabilization of intermedi- 
ate compositions for both the solid and the liquid, and 
so a non-linear interpolation would be required to detect 
enhanced thermostability from the background mixing 
effects. 

In order to understand better how particular struc- 
tures are stabilized, it is useful to decompose the energy 
into different terms. Firstly, we divide the energies into 
contributions from those pairs of atoms that are nearest 
neighbours from those that are not, i.e. 
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where E^n and Ennn are the energetic contributions from 
nearest neighbours and non-nearest neighbours. Nearest 
neighbours are simply defined using a distance criterion. 
In practice we used 1.3 times the relevant equilibrium 
pair separation, but as there is usually a clear separation 
between nearest-neighbour and next-neighbour shells for 
the ordered structures that we find to be the global min- 
ima, the precise value is not critical. 

Secondly, we divide up E^n by defining a strain energy 
in terms of the difference in energy between the nearest- 
neighbour energy, and the energy this term would have if 
all the nearest neighbours had the same separation (mea- 
sured with respect to the relevant equilibrium pair sepa- 
ration): i.e. 
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where n^n is the number of nearest neighbours, and r'^,-,, 
the reduced average separation of these nearest neigh- 
bours is defined as {rij/aap)nn, where the average is 
taken over only those pairs of atoms that are near- 
est neighboursi^ Hence, ii^strain measures the energetic 
penalty arising because a structure has a distribution of 
nearest-neighbour distances. 

As -Bnnn IS relatively insensitive to structure, the 
best structure is determined mainly by the balance be- 
tween maximizing nnn and minimizing i^strain- For one- 
component polytetrahedral clusters, both rinn and i^strain 
are usually large, and by following -Estrain as a function 
of (Jbb/'^aa will allow us to see how this strain can be 
relieved by having different sized atoms. 



B. Global optimization 

The basin-hopping global optimization algorithm^i^ 
has proved to be very successful for the optimization of 
a wide range of clusters. The method simply involves 



doing a constant temperature Metropolis Monte Carlo 
simulation, where after each step a local minimization 
is performed on the resulting configuration, and the ac- 
ceptance criterion is based upon the energies of these 
minimized configurations. The success of the method 
is because it effectively searches a transformed potential 
energy surface, for which the dynamics and thermody- 
namics is more helpful for optimization^2iiS 

In addition to the standard basin-hopping moves (ran- 
dom displacements of all atoms and rotations of low- 
energy atoms around the centre of mass) used for one- 
component systems, for binary clusters it is important 
to incorporate moves into the algorithm that allow effi- 
cient exploration of the space of homotops and compo- 
sition. The two moves that we used involved swapping 
the identities of an A and a B atom^iSi and changing the 
identity of a single atom. In optimizations runs where the 
composition was held fixed, only the former composition- 
preserving moves were used. Typically, these types of 
moves represented 50% or more of the total moves. 

The current optimization task of finding the structure 
and composition of a binary cluster that has the lowest 
energy (we shall call this the 'compositional global min- 
imum') is very challenging, because the number of pos- 
sible structures increases very rapidly with size. Firstly, 
the number of geometric isomers scale as exp(aA'^)i^ii^ 
Secondly, the number of possible homotops for a given ge- 
ometric structure is proportional to (jy ) (assuming that 
all the homotops are locally stable). Thirdly, the num- 
ber of possible compositions increases linearly with TV. 
For example, for LJioo the estimated number of geomet- 
ric isomers is of the order of 10'^^. If all these are stable 
for the binary clusters and all homotops are possible, this 
then gives an estimate of 10^^ energetically different min- 
ima in the space that we search to find the compositional 
global minimum of BLJioq. It is, therefore, important 
that we make as much use of the structural information 
gleaned from unbiased runs to aid the optimization task. 

Our initial strategy was to perform an initial scries of 
runs for each cluster size from random starting points at a 
moderate temperature, similar to that typically used for 
one-component LJ clusters. Subsequently, we performed 
low temperature runs from the best structures found from 
the above runs. The rationale is that for gbb/'^aa values 
sufficiently close to one, the energy differences between 
homotops can be significantly smaller than that between 
geometric structures. For example, this feature mani- 
fests itself in the thermodynamics of such binary clus- 
ters as low-temperature heat capacity peaks associated 
with permutational disordering. '^'^■'*'^ The hope is then 
that the high temperature runs would find the best geo- 
metric structures, and the low temperature run the best 
homotop consistent with this geometric structure. At 
larger obb/'^aAi this separation of energy scales breaks 
down, and so the low-temperature runs are less impor- 
tant. We also performed series of runs, where the best 
structure for a particular size and size ratio were used as 
starting points for basin- hopping runs at different size ra- 



tios and sizes (with the requisite number of atoms added 
or removed) . These different types of runs were appHed 
iteratively until convergence seemed to be reached. 



III. RESULTS 



Case studies 



Before we present results whilst directly searching for 
the global minima in both composition and configura- 
tion space, it is useful to get a feel for the energetics 
as a function of composition. To do this, we map out 
the landscape associated with the energy of the global 
minima as a function of ha and ubb I o'aa for a number 
of examples, namely BLJ13, BLJ45 and BLJ55. For the 
latter two, this represents a significant effort, since we 
have found the global minima at 1320 and 1620 different 
points in this space, respectively. 

For BLJ13, the structure of the global minima is al- 
ways the centred icosahedron, so Fig. |2Ia) simply rep- 
resents the variation of the energy of this one structure. 
Close to obbI^aa = 1, as expected, the structure with 
one small atom at the centre of the cluster is the most 
stable composition. The distance between adjacent ver- 
tices of the perfect icosahedra is 5.15% longer than that 
between a vertex and the centre. Consistent with this, 
the strain energy goes to zero when obbI <y ab takes this 
value, i.e. obbI <^ aa — 1.1085. The minimum in the total 
energy is displaced to slightly larger obbJoaa, namely 
1.1303, because this leads to a greater contribution to the 
energy from next neighbours. Beyond this, the energy 
of the AB12 structure rises, and at aBB/<^AA = 1.1614, 
the A4B9 structure becomes the most stable composition. 
With three smaller atoms in the surface, this structure is 
more able to relieve the compressive strain that builds 
up in the AB12 structure. With further increases in 
^^bbI <^AA the global minimum changes twice more with 
the number of A atoms increasing in order to achieve a 
lower strain energy. These transitions give a flat bottom 
to the landscape in Fig. |2a). 

BLJ45 shows a much richer structural behaviour, be- 
cause there are a variety of competing geometric struc- 
tures. Firstly, there are two ways that atoms can be 
added around the 13-atom Mackay icosahedron, as illus- 
trated in Fig. 121 In general, the Mackay overlayer con- 
tinues the face-centred-cubic (fee) packing of the twenty 
vertex-sharing fee tetrahedra that make up a Mackay 
icosahedron, whereas the anti-Mackay overlayer adds 
atoms in sites that are hep with respect to the fee tetra- 
hedra. The anti-Mackay layer has a lower surface density, 
and consequently has a greater number of nearest neigh- 
bours, but for monoatomic systems also a larger strain 
energy. Typically, therefore, growth starts off in the anti- 
Mackay overlayer, and then switches to the Mackay over- 
layer, as the layer grows. 

For growth on the 13-atom icosahedron, the anti- 
Mackay overlayer maintains the polytetrahedral charac- 
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FIG. 2: (Colour online) BLJ13. (a) Landscape showing the 
dependence of the energy of the global minimum on (tsb/ctaa 
and Na- (b) The four compositional global minima and the 
dependence of (c) their total energies and (c) strain energies 

on CTBB/f^AA- 



ter of the structures with each interior atom having a 
local icosahedral coordination shell (we call such struc- 
tures polyicosahedral) , whereas the Mackay overlayer in- 
troduces some octahedral interstices into the structure. 
For LJ clusters, A^ = 31 is the first size at which the 
global minimum adopts the Mackay overlayer.^ We chose 
to examine BLJ45 in detail, because at this size the anti- 
Mackay overlayer is able to be completed, giving a struc- 
ture with point group Ih that we call the anti-Mackay 



Anti-Mackay 
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FIG. 3: The two possible overlayers for growth around a com- 
plete Mackay icosahedron. These are illustrated for a single 
face of the underlying 13-atoni (top) and 55-atoni (bottom) 
Mackay icosahedron. 



icosahedron. 

Secondly, there is the possibility of structures that are 
polytetrahedral, but that are not polyicosahedral. In- 
stead, there are interior atoms that have a coordina- 
tion number (Z) greater than 12. For polyicosahedral 
structures, there are five tetrahedra around every inte- 
rior edge. However, as we noted in the introduction, the 
corresponding packings of regular tetrahedra involve gaps 
that increase rapidly with the size of the packing (Fig. 
^. As a consequence, it is not possible to generate a 
bulk polyicosahedral packing, and for clusters, the sur- 
face density decreases and the average nearest-neighbour 
distance between the surface atoms increases as the struc- 
tures become bigger. 

An alternative way to pack tetrahedra is to intro- 
duce some edges in the structure that have six tetrahe- 
dra around them. Such an arrangement is locally less 
favourable — for regular tetrahedra there is an overlap 
equivalent to 63.17°. However, for sufficiently large pack- 
ings, the combination of 'gaps' and 'overlaps' allows one 
to generate structures that have less overall strain than 
polyicosahedral packings of the same size. For exam- 
ple, bulk polytetrahedral structures, called Frank-Kasper 
phases ;^ii^2iH are now possible with these two types of 
edge environments. For clusters, such polytetrahedral 
structures have a higher surface density than polyicosa- 
hedral structures, and so need a smaller atomic size ratio 
to be stabilized in binary systems. 

Those edges that have six tetrahedra around them are 
said to have disclination lines running along them, and 
such polytetrahedral structures can be can be viewed 
in terms of a network of disclination lines threading an 
icosahedrally-coordinated medium. Those atoms that 
have disclinations passing through them have a coordina- 
tion number greater than 12. For example, 14-coordinate 
atoms have a single disclination line passing through 
them, and 15- and 16-coordinate atoms act as nodes for 
three and four disclination lines, respectively. We will 



use the presence of these different types of Frank-Kasper 
coordination polyhedra to differentiate these polytetra- 
hedral structures. 

The [aBB/^AA^NA) landscape associated with the 
BL J45 global minimum is shown in Fig. 0] It is notewor- 
thy that the there is no optimal size ratio, but instead the 
energy of the global minimum decreases virtually mono- 
tonically with increasing size disparity. Also apparent is 
the magnitude of the stabilizations that can be achieved 
compared to monoatomic LJ clusters. The lowest-energy 
structure in the ctbb /(^aa range that we consider here is 
27.252e lower in energy than the LJ45 global minimum, 
i.e. a 12.7% decrease in the total energy. 

There are four different compositional BLJ45 global 
minima in the parameter range studied. Near to 
(^bb/'^aa = 1, as one would expect the global minimum 
is icosahedral with a Mackay overlayer, and has the same 
geometric structure as for LJ clusters. The ideal compo- 
sition of this structure depends slightly on <tbb /(^aa- At 
(^bb/'^aa — 1.057 the polytetrahedral structure AgBsy 
with a single negative disclination threading the struc- 
ture becomes most stable. At ubb/i^aa = 1.071 the 
core-shell anti-Mackay icosahedron becomes most stable. 

These changes are driven by a balance between the 
number of nearest neighbours and the strain energy of 
the structures. As (Jbb/'^aa increases, structures with 
a larger number of nearest neighbours become the most 
stable, when their strain energy is sufficiently reduced. 
In this series, the number of nearest neighbours increases 
from 180 to 192 to 204, but at small obbI'^aa so does the 
strain energy also increase. In particular, notice that the 
strain energy for the polytetrahedral structure involving 
disclinations is less than the polyicosahedral structure 
(Fig.Hb)). 

These results illustrate how allowing two atoms of dif- 
ferent sizes can substantially reduce the strain energy 
associated with polytetrahedral structures. For exam- 
ple, the strain energy of the anti-Mackay icosahedron 
can be decreased from 24.4 to 1.7e (Fig^jj). However, 
unlike BLJ13, the strain cannot be completely removed 
There are only two different nearest-neighbour distances 
in a 13-atom icosahedron, whereas there are six different 
nearest-neighbour distances for the anti-Mackay icosahe- 
dron, and they cannot be adjusted so that they all have 
the same Tij/aap- 

When we also look at the compositional dependence 
of the global minimum in Fig. El^b), we should remem- 
ber that at Na = 0, Na — N and ctbb/^aa — 1 the 
model is the same as the one-component Lennard- Jones 
model. Similarly, close to these values the structure has a 
Mackay overlayer, as for LJ clusters, i.e. for clusters with 
mostly A, mostly B or small <jbb I (^aa- By contrast, the 
polyicosahedral structures are most stable for intermedi- 
ate compositions and larger gbbI'^aa- For much of the 
structural phase diagram, sandwiched between these two 
structural types is a region where polytetrahedral struc- 
tures with disclinations are most stable. However, there 
is some asymmetry, as at larger ubbI'^aa-, this zone is 
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FIG. 4; (Colour online) BLJ45. (a) Landscape showing the dependence of the energy of the global minimum on obbJoaa and 
Na- (b) Zero-temperature structural phase diagram showing how the structure of the global minimum depends on ubb /o^aa 
and Na, superimposed on the contour plot of the landscape in (a). Each data point corresponds to a change in the structure of 
the global minimum as a function of (tbb/cfaa- The region where the anti-Mackay icosahedron is global minimum is labelled Ih- 
(c) The four compositional global minima and the dependence of (d) their total energies and (e) strain energies on (jbb /uaa- 



only found for larger values of Na- 

Our final case study is BLJ55. This size has been 
chosen, as it corresponds to the size at which a com- 
plete Mackay icosahedron is possible. Again, we illus- 
trate the landscape associated with the global minimum 
and, in this case, the seven compositional global min- 
ima (Fig. [SJ. The overall behaviour is somewhat simi- 
lar to BLJ45. There are four main differences. Firstly, 
the Mackay structures are most stable over a greater 
range of (JbbI<^aa and composition, which is unsurpris- 
ing given that ./V = 55 is a Mackay magic number. Sec- 
ondly, there is a more pronounced change in slope in 
Fig. Efa) associated with the onset of polytetrahedral 
structures. The landscape is relatively flat in the re- 
gions where Mackay structures are most favoured, but 
goes down more steeply, in the region where polytetra- 
hedral structures are most stable, because the strain en- 
ergy associated with these structures rapidly decreases, 
as the size disparity increases (Fig.|SJe)). Thirdly, polyi- 
cosahedral structures are most stable for a smaller range 
of parameters. This is because, as the size increase the 



strain energy associated with these structures for the 
monoatomic LJ clusters increases rapidly, and hence the 
size difference that is needed to stabilize them also in- 
creases. Fourthly, as a corollary of the above, the area 
in Fig. Ela) where polytetrahedral structures involving 
disclinations are more stable increases. Furthermore, 
there are two types of such structure, which can be differ- 
entiated by the nodes in the disclination network. The 
'Z14' structures involve a single disclination line pass- 
ing through the cluster, whereas the 'Z15' structures has 
three disclinations radiating out from the 15-coordinate 
atom. 

For BLJ55, the compositional global minimum changes 
six times as <tbb / c^aa increases. First, there is the core- 
shell Mackay icosahedron. Although the strain in this 
structure is much less than for polytetrahedral structures, 
the strain energy can still be significantly reduced com- 
pared to the LJ case by introduction of different-sized 
atoms (Fig. IHe)). For example, this structure has its 
lowest energy at obbI <^ aa = 1-064. The next three com- 
positional global minimum are polytetrahedral structures 
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FIG. 5: (Colour online) BLJ55. (a) Landscape showing the dependence of the energy of the global minimum on obbJoaa and 
Na- (b) Zero-temperature structural phase diagram showing how the structure of the global minimum depends on ubb /o^aa 
and Na, superimposed on the contour plot of the landscape in (a). Each data point corresponds to a change in the structure 
of the global minimum as a function of (Jbb/o-aa- The region where polytetrahedral structures involving disclinations are the 
global minima is divided into two subregions, depending on whether the highest coordination number for an atom is 14 or 15. 
(c) The seven compositional global minima and the dependence of (d) their total energies and (e) strain energies on obb/ctaa- 
In (c) for the relevant global minima the disclination networks are depicted to the left of the structure. 



that involve disclination lines. Note, that the Z\^ struc- 
ture is more stable at smaller ubbI'^aa as it involves 
a greater disclination density than the Z14 structures. 
These structures are not pure core-shell clusters, as al- 
though all the 12-coordinate interior atoms are A atoms, 
there is a preference for the larger atoms to lie on the 
disclinations, which is unsurprising given their higher co- 
ordination number. The final three compositional global 
minima are polyicosahedral, and differ just in the posi- 
tion of capping atoms and composition. 



B. Compositional global minima 

In this section, we focus our attention on the compo- 
sitional global minima. Figs. El [7| |H1 and|51 summarize 
the results. The energies and points files for all the pu- 



tative global minima are available online44 In Fig. El the 
energies at different size ratio are compared. This figure 
clearly illustrates that, as noted in the above case stud- 
ies, the energy virtually monotonically decreases with in- 
creasing size ratio. It also shows the substantial nature 
of the stabilizations that are achieved compared to the 
one-component Lennard-Jones system. 

In Fig.[7| the energies at the six size ratios we consider 
are plotted in a way that reveals the particularly stable 
sizes. In Fig. |Slwe show some of the particularly inter- 
esting or stable structures associated with the different 
structural types. 

Fig. provides a structural phase diagram showing 
how the (N ^a BB I <^ aa) plane can be divided into regions 
where the compositional global minima have the same 
type of structure. To construct this diagram, we consid- 
ered intervals of 0.01 in ctbb/'^aa and reoptimized the 
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FIG. 6: The energy of the BLJ global minimum for the six val- 
ues of (JbbIoaa studied, relative to E^^i ^ fit to the energies 
of the Mackay icosahedra for LJ clusters. A line correspond- 
ing to the LJ global minima is also included 



five lowest-energy structures from the nearest size ratios 
for which global optimization was performed. Then, we 
checked if the structure of the global minimum changes 
in any of these intervals, and if so we obtained the pre- 
cise value of obbJoaa at which this change took place. 
Of course, this approach will potentially miss the true 
global minima at intermediate values of ubb I<^aa if it is 
not one of the five best for the values at which we ran 
the global optimization algorithm, and it only allows for 
one change in the global minimum in any 0.01 interval. 
However, we are concerned more with the overall form 
of this diagram, on which these approximations will only 
have a very minor effect. 

The reference system to which to compare our re- 
sults is, of course, the one-component LJ clusters, for 
which structures based on the Mackay icosahedra are 
dominant in the current size range. Growth around the 
13- atom icosahedron initially occurs in the anti-Mackay 
overlayer, but for LJ31 and beyond the global minimum 
has a Mackay overlayer.— Similarly, for the growth of the 
next icosahedral shell, the global minima initially have 
an anti-Mackay overlayer, but for LJ82 and beyond (with 
the exception of LJgs) the Mackay overlayer is more sta- 
ble. The only exceptions to this dominance of icosahe- 
dral structures are for LJ3g, LJ75„77 and LJgg for which 
an fee truncated octahedron, Marks decahedra'''^ and a 
Leary tetrahedron^ respectively, are just more stable 
than the competing icosahedral structures. 

As one moves away from gbb /<^aa = 1 in the struc- 
tural phase diagram, firstly icosahedral structures quickly 
become more stable than the five non-icosahedral struc- 
tures mentioned above. Secondly, for both the second 
and third icosahedral shells the crossover size at which 
the Mackay overlayer becomes more stable is pushed to 
larger sizes (Fig. O. This is because the anti-Mackay 
overlayer is more strained, and so the introduction of two 
atomic sizes can lead to a greater reduction in its strain 
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FIG. 7: The seven panels correspond to the energies of the 
global minimum for the six values of (tbbJoaa studied and for 
the one-component LJ clusters, relative to Ea,-je{aBB/(yAA), a 
fit to the energies of the global minima at that size ratio using 
the form a + hN^'''^ + ciV^/^ + dN . 



energy. For both types of overlayer, the larger atoms will 
preferentially go into the surface layer, because of the ten- 
sile strain in the surface of the icosahedral LJ clusters, 
leading to the formation of core-shell clusters, as already 
illustrated by A13B32 and A13B42 in the previous sec- 
tion. There are many further examples in Fig. |H1 As for 
the 13-atom icosahedron considered in the last section, 
there is an optimal size ratio for these core-shell struc- 
tures. Beyond this, compressive strains begin to build up 
in the surface until at some point it becomes favourable 
to include some of the smaller A atoms in the surface. 




FIG. 8: (Colour online) A selection of the particularly stable BLJ global minima grouped according to their structural type: (a) 
polyicosahedral, (b) polytetrahedral with disclinations and (c) the 55-atom Mackay icosahedron with an anti-Mackay overlayer. 
For the larger clusters, to the right of the cluster, the A-atom core is also depicted, and in (b) to the left is the disclination 
network. 



Some examples of the structures that result are shown 
in Fig. IHta), e.g. A7B12, AioBig and Ai2B22- However, 
the window of ctbbI'^aa values for which these core-shell 
structures are most stable is wide, illustrating their par- 
ticular stability. 

This preference for the anti-Mackay overlayer can also 
be seen in the observed magic numbers (Fig.[7J). For the 
second shell, anti-Mackay clusters are particularly stable 
at Af=19, 23, 26 and 29 for LJ clusters, and correspond to 
2, 3, 4 and 5 interpenetrating icosahedra. These become 
more prominent as obbJoaa increases with additional 
magic numbers at 7V=34 and 45, and to a lesser extent 
at iV=32, 37 and 39. The magic number at TV = 45 
is particularly strong as it corresponds to the comple- 
tion of the anti-Mackay overlayer. Interestingly, beyond 
A'' = 34 the way the anti-Mackay overlayer grows de- 
pends on OBBl<yAA- At smaller ubb/o'aa the formation 



of structures that can be thought of as interpenetrating 
complete icosahedra continues, as illustrated by A9B30. 
In contrast, at larger (Jbb /(^aa the sites above the faces 
of the central 13-atom icosahedra are completely filled 
and further growth occurs just by adding atoms above 
the vertices (Fig. |3J), e.g. A13B26. 

Similarly, for the third shell, as well as the magic num- 
ber at N=71 already present for LJ clusters, additional 
magic numbers appear at 81, 88 and 95. These corre- 
spond to covering five, eight, ten and twelve faces of the 
underlying Mackay icosahedron, respectively. 

As already mentioned, clusters with an anti-Mackay 
overlayer covering the 13-atom icosahedron are polyicosa- 
hedral in character. However, there is no reason why this 
type of packing should not continue beyond the comple- 
tion of this overlayer at A^ = 45. Indeed, for the range of 
o'bb/o'aa values considered here such structures are pos- 
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the Z15Z16 structures, e.g. A26B61, have one node where 
three discUnations meet, and one where four meet. 

From Fig. O some systematic trends in the charac- 
ter of the dischnation networks are clear. As N in- 
creases, the structures next to the Mackay icosahedral- 
polytetrahedral boundary have an increasing number of 
dischnations in order to reduce the growing strains that 
would have otherwise occurred. However, as (Tss/cyiA 
increases the different-sized atoms are able to reduce 
some of this strain, and the disclination density goes 
down, until the boundary with disclination-free polyi- 
cosahedral structures is reached. So, for example, for 
BLJ82 the structure changes from 2Z16, to Z16Z15 to 
2Z15 to Zlb to ZlA as (Tbb/<^aa increases. 



FIG. 9: (Colour online) Structural phase diagram showing 
how the structure of the global minimum depends on A'' and 
obb/i^aa- Each point corresponds to the value of (jbb/o-aa 
at which the global minimum for a given N changes. The lines 
divide the diagram into regions where the global minima have 
the same structural type. The labels 7i3 and Jss stand for the 
13- and 55-atom Mackay icosahedra. 



sible up until N = 69, and the magic numbers at N—bQ^ 
62 and 66 that appear at larger aBBl<^AA correspond 
to core-shell polyicosahedral structures with the double, 
triple and quadruple icosahedra as their core (Fig.EJa)). 
Also interesting is the A23B41 structure, for which the 
core of A atoms is two face-sharing icosahedra. 

In the structural phase diagram as N increases a re- 
gion opens up between the Mackay icosahedral and poly- 
icosahedral structures, where polytetrahedral structures 
involving dischnations are most stable. The first such 
structure occurs at A^ = 44 and as N increases, these 
structures occupying an increasing proportion of the 
phase diagram. These structures are not seen at small 
size, firstly because the strain associated with polyicosa- 
hedral structures is not prohibitively high, and secondly 
because the introduction of dischnations would lead to 
high disclination densities, and unfavourable strain en- 
ergies. For this reason the Frank-Kasper coordination 
polyhedra that are possible at N—\h^ 16 and 17 are never 
most stable. 

Later when substantial strains have builit up in the 
polyicosahedral structures, the introduction of dischna- 
tions along a small minority of the edges can lead to 
a reduction in the overall strain energy. The struc- 
tures A12B39, A12B42 and A12B49 provide interesting 
examples, and correspond to the covering of the above 
Frank-Kasper coordination polyhedra by a complete anti- 
Mackay-like overlayer with atoms added above every face 
and vertex (except the six- fold vertices for the 51- and 
54-atom structures). 

At larger sizes more complex disclination networks are 
possible. For example the 2Z16 structures, e.g. AigB64, 
A30B62, have two nodes where four dischnations meet, 
and gives rise to an ethane-like disclination network, and 



IV. CONCLUSIONS 

In this paper we have introduced a new global opti- 
mization approach for binary clusters that for a given 
size searches directly for the composition of greatest sta- 
bility. This development makes it tractable to explore the 
size evolution of binary cluster structure systematically, 
rather than just arbitrarily selecting particular sizes or 
compositions. Of course, by focussing on the composi- 
tional global minimum, some information is lost, but as 
one is usually most interested in identifying the most sta- 
ble magic number clusters, this is of little import. Besides 
the alternative approach of finding the global mininum 
at each size and composition would both be computa- 
tionally extremely challenging and result in a surfeit of 
information, the added value of which is far from clear. 
The approach is straightforward to apply to other binary 
systems, and, for example, we have just completed ap- 
plying this approach to all clusters with up to 150 atom 
for a number of binary metal systems*^ 

One of the most important aims of theoretical studies 
of structure is to provide models that can aid the inter- 
pretation of experimental observations. In this regard, it 
is reasonable to ask whether our approach of optimizing 
the composition is realistic of what might be occurring 
experimentally. The answer would of course depend on 
how the binary clusters are produced experimentally, but 
recent experiments on copper-tin clusters showed that 
the particularly stable compositions could be obtained 
after annealing the clustersi^ 

We have used the current approach to systematically 
explore the structural effects of having different-sized 
atoms on binary LJ clusters. Particularly interesting is 
how this stabilizes polytetrahedral clusters, both those 
that are polyicosahedral and those that involve disch- 
nations, which for the one-component LJ clusters are 
too strained to be competitive. The effects of the two 
atom sizes is somewhat similar to other methods of strain 
relief, such as widening the potential well, which also 
pushes the anti-Mackay to Mackay transitions to larger 
size and stabilizes polytetrahedral structures with discli- 
nation lines ii^iii"?i In future work, we plan to extend our 
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exploration of the BLJ system to other choices for the 
four parameters in the potential, and to help develop its 
status as an archetypal model system for which to un- 
derstand the structure of binary clusters. 

Our results have provided a zoo of interesting struc- 
tures, which are good candidates for particularly stable 
clusters for binary systems where the two atom types 
have significant differences in size. Although BLJ clus- 
ters are not a realistic model for much, except for rare gas 
mixtures at certain parameter choices, as for monoatomic 
LJ clusters, our expectation is that the stable structural 
forms seen for the BLJ clusters are robust, and likely to 
be relevant for a wide variety of systems, where the in- 
teractions are approximately isotropic. This confidence 
is based on previous experience, where particular struc- 
tural forms first seen for a model system are later found 
in real systems. For example, the Leary tetrahedrouj^ 
which was first discovered as the global minimum for 
LJgg, was later found to be one of the dominant struc- 
tural forms for clusters of Ceo moleculesi^ Indeed, core- 
shell polyicosahedral structures have already been found 
for binary metal clusters, such as the Ag-Ni and Ag-Cu 
systems .Mii^ilL^ 

In this paper, we have focussed on the global minima 
of our system, and so Fig. represents the structural 
phase diagram at zero temperature. It is, therefore legit- 
imate to ask how our results would be affected by tem- 
perature. Firstly, for sufficiently small size differences 
between the atoms, it is likely that there will be low- 
temperature order-disorder transitions, where the geo- 
metric structure is retained, but the permutational order 
is losti^i^ However, as gbbJc^aa increases, these tran- 
sitions are likely to become less common, because the A 
and B atoms quickly develop strongly preferred positions, 
as, for example, in the core-shell ordering. 

Secondly, it is known for the L J system that an effect of 
temperature is to push the anti-Mackay to Mackay transi- 
tion to larger size^^SiSi^iiS because of the greater vibra- 
tional entropy of the anti-Mackay structuresiii Indeed, 
the form of Fig. |2| at small (JbbI<^aa looks similar to 
the [N, T) structural phase diagram for LJ clustersiSi^ 
Therefore, temperature will reinforce the stabilization of 
polytetrahedral structures seen for BLJ clusters. 

Here, we have probed all BLJ clusters with up to 100 
atoms. But what would one expect at larger sizes, and 
in the bulk limit? The dominance of polytetrahedral 
structures in the current size range is because of their 
favourable surface energetics, and it can be seen in Fig. 
Elthat the range of stability for non-polytetrahedral clus- 
ters is increasing as N increases. The latter would sug- 
gest that Frank-Kasper phases are never the ground state 
in the bulk limit. Instead, phase-separated A and B fee 
crystals are likely to have the lowest energy for bulk. 



However, one cannot rule out that this boundary might 
flatten off at larger sizes, nor that there might be some 
range of temperature and pressure, where Frank-Kasper 
phases might be most stable. 

The reverse of the above prediction is that there are 
likely to be binary systems, which although they do not 
have any stable bulk Frank-Kasper phases, neverthe- 
less exhibit Frank-Kasper-like structures for their clus- 
ters. Indeed, we have already identified such struc- 
tures for Ni-Al, Ag-Cu and Ag-Ni clusters using the 
current methodology^ Although likely to be less com- 
mon than for binary clusters, there are a few one- 
component systems that also exhibit these polytetrahe- 
dral clustersi2iSS*S 

Our results are also of relevance for understanding the 
effect of size mismatch on glass formation. In order to 
explain the difficulty of homogeneous nucleation in super- 
cooled metal droplets, Frank argued that if the preferred 
local order within the liquid was incompatible with the 
local structure in the crystal, nucleation would be sup- 
pressed because of the substantial structural rearrange- 
ment requiredi^ Furthermore, he used the structures of 
isolated LJ clusters, in particular the stability of the 13- 
atom icosahedron compared to an fee cluster of the same 
size, to provide a picture of the local order within the 
liquid. More recently, the predictive power of isolated 
clusters for understanding liquid structure has received 
significant empirical support i^^iffiiSi 

Empirically, it has been found that one of the con- 
ditions for metallic alloys to form bulk glasses is that 
there is a size difference of at least 12%iSS A simple the- 
oretical justification of this can be given in terms of the 
size difference required to destabilize a crystalline solid 
solutionSSi For the current system, the critical value of 
(^bb/^^aa at which a solid solution is no longer always 
stable has been found to be 1.20.^''^ Although the deter- 
minants of glass-forming ability are subtle^SLSi and more 
sophisticated theories can be developed^SiSSiZS our re- 
sults highlight the potential role played by the local struc- 
ture within the liquidj^i and how size differences can sta- 
bilize local polytetrahedral, in particular icosahedral)^^ 
order, and hence frustrate crystallization. Indeed, one 
of the commonly used model glass-forming systems used 
in the simulation community is a 50:50 BLJ mixture 
with ubb/'^aa — 1-2, and all other parameters as in 
the present studyiSiZiZk 
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